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ABSTRACT 

Supermassive black hole binaries (SMBHBs) are the products of frequent 
galaxy mergers. The coalescence of the SMBHBs is a distinct source of gravita- 
tional wave (GW) radiation. The detections of the strong GW radiation and their 
possible electromagnetic counterparts are essential. Numerical relativity suggests 
that the post-merger supermassive black hole (SMBH) gets a kick velocity up 
to 4000 km s -1 due to the anisotropic GW radiations. Here we investigate the 
dynamical co-evolution and interaction of the recoiling SMBHs and their galac- 
tic stellar environments with one million direct iV-body simulations including 
the stellar tidal disruption by the recoiling SMBHs. Our results show that the 
accretion of disrupted stars does not significantly affect the SMBH dynamical 
evolution. We investigate the stellar tidal disruption rates as a function of the 
dynamical evolution of oscillating SMBHs in the galactic nuclei. Our simulations 
show that most of stellar tidal disruptions are contributed by the unbound stars 
and occur when the oscillating SMBHs pass through the galactic center. The av- 
eraged disruption rate is ~ 10~ 6 M Q yr -1 , which is about an order of magnitude 
lower than that by a stationary SMBH at similar galactic nuclei. Our results also 
show that a bound star cluster is around the oscillating SMBH of about ~ 0.7% 
the black hole mass. In addition, we discover a massive cloud of unbound stars 
following the oscillating SMBH. We also investigate the dependence of the results 
on the SMBH masses and density slopes of the galactic nuclei. 
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INTRODUCTION 



Supermassive black hole binaries (SMBHBs) are pre dicted by the hierarchical galaxy for- 
matio n model in A cold dark matter (ACDM) cosmology ( Begelman et alll980 : Volonteri et al 
20031 ). For merging galaxies, their two SMBHs with galactic cores will firstly approach each 
other by dynamical friction, and then get close enough to form a bound, compact binary 
system. After that, the SMBHB ma y stall at a so called "h ard binary separation" for a time 
even longer than the Hubble time (IBegelman et al.lll980l ). However, recent investigations 
suggest that the hardening rates of SMBHBs can be boosted and they may coalesce within 
the Hubble time either due to various stellar dynamical processes other than spherical two 



body relaxation ([Yull2002l : IChatterjee et al 



2003 



Merritt fc Poon 



2004 



Sesana et al. 20081 ) . or gas dynamics ( Gould fc Rix! 200ol : IColpi fc Dotti 
therein) . 



Berczik et al. 2006 



2009 



and references 



Since the evolution of SMBHBs deeply impacts the evolution of host galaxies, it is very 
important for us to find observational evidences to constrain evolution models of SMBHBs. 
A statistic way is the calculation of tidal disruption rates. A dormant SMBH could be 
temporarily a ctivated by tidally disrupting a star passing by a nd accreting the d isrupted 



1989 



stellar debris (IHillsl Il975t iReesl Il988t lEvans fc Kochanek 
may have been observed in se ve ral non-acti v e gala xies ( Komossa Sz Bade 



Lodato et all l2009h. which 



20021 : iGezari et all 120081 120091 ). IChen et al.l (120081 ) and IChen et all (120091 ) calculated the 



1999 



Komossa 



tidal disruption rate in SMBHB systems at different evolutionary stages, and found that is 
significantly different from the typical rate in a single SMBH for several orders of magnitude. 

The most certain proof for detecting SMBHB individually may come from gravitational 
wave (G W) radiation observation. Coalesc ing SMBHBs are important sources of GW ra- 
diation (jPetersI Il964t IBegelman et al.lll980[ ). and can be detected within the next decades 
by the L aser Interferometer Space Antenna (LISA) and the Pulsar Timing Array (PTA) 



program (IBerentzen et al.ll2009 



Amaro-Seoane et al.ll2010l ). Because of the poor spatial res- 



olution of both LISA and PTA for locating GW radiation sources, it is of key importance to 
detect electromagnetic counterparts (EMCs) of GW radiation sources. Besides, identifying 
SMBHBs by their EMCs is also essential to constrain the poorly understood galaxy-merger 
history. 



Several EMCs have been suggested in the literature to probe SMBHBs, (1) preces- 
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sion of jet orientation and i ts acceleration in radio galaxies during the in-spiraling of SMB- 
HBs f IBegelman et al.lll980l ; iLiu fc Chenl 120071 ). (2) o ptical periodic outbursts in blazars due 



to th e interaction of S MBHBs and accretion disk (jSillanpaa et al.l Il988l ; ILiu et al.l 11995 



20061 : ILiu fc Wu! 120021 : IValtonen et al.l 120081 : lHaiman et al.l 120091 ). (3) jet reorientation in 



X-shaped radio gal axies due to the exchange of angular momentum between SMBHBs 
and a ccretion d i sk (|Liul 120041 ). (4) interruption of tidal disruption flares in SMBHBs sys- 
tems (ILiu et al.l 120091 ). Besides, there are also some EMCs to probe the coalescence of 

intermittent activity in double-double radio galaxies at 
20"03l). R afterglow following; 



SMBHBs and its remnant, 



binary coalescence ( Liu et al 



f lMadau fc Quataertl 12004 



binary coalescence (iMilosaylievic fc Phinnevl l2005t IShields fc Bonningj l2008t iLippai et al. 



20081 : ISchnittman fc Krolikl 12008). (3) system atically shifted broad emission lines relative 



to narrow emission lines (Komossa et al.ll2008l ) and off-center active galactic nuclei (AGNs 



Loebll2007n because of SMBH GW radiatio n recoil, (4)tidal dis- 



ruption flares and hypercompact stellar systems for recoiling black hol es (IKomossa &: Merritt 
20081 : lO'Learv fc LoerJl20~0a IMerritt et alibood : Istone fc Loeblboilh . 



The breakthrough on numerical rela t ivity in the past few years reveals the fina l stage 
of SMBHB's coalescence JPretoriusI [2OO5I : ICampanelli et al.ll2006[ baker et al.ll2006ah . The 
coalescence remnant SMBH can be recoiled due to the anisotropically GW emission during 
the inspiral and fin al coalescence JPeres 1962 ). For nonspinning SMBHs, the recoil ve locity is 
Vkkk < 200 km s" 1 faaker et al.l l2006bl : [Gonzalez et al.1l2007al : Irlerrmann et al1l2007f ). which 
is just as the same order of magnitude as the stellar dispersion velocity in the galactic center. 
However, if both of two SMBHs are rapidly spinning, the recoil velocity can be as large as 
thousands kilometers per second. This recoil velocity depends sensitively on the intersection 
angle between spin vectors of two SMBHs and their linear momenta. For some extreme 
cases, with maximally spinning equal-mass SMBHs and antialigned spins in orbital plane, 
the re coiling velocity can achieve to ~ 4000kms _1 (ICampanelli et al.l 120071 ; I Gonzalez et al. 
2007bh . 



A recoiling SMBH with high velocity has significant displacement re lative to the galacti c 
center. For most of massive galaxies with escape velocity < 3000 kms" 1 (IMerritt et al.ll2004l ). 
the recoiling SMBH has opportunity to escape from its host galaxy. Those recoiling SMBHs 
which are still bound to host galaxies will oscillate around the galact ic centers and change 
the core density profiles of the ho st galaxies due to dynamical friction ( iBoylan-Kolchin et al. 
20041 : iGualandris &; Merrittll2008l ). In addition, the recoiling SMBH with a fr action of its ac 
cretio n disk can produce sets of emission lines sep arated by relative velocity (IKomossa et al. 
20081 ) or off-center active galactic nuclei (AGNs) ( IMadau fc Quataertl 12004 ; iLoeblboOTT ) . 



In gas poor environment, a recoiling SMBH might be observed as a "hypercompact stel- 
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lar system" (HCSS). That is a kind of star cluster which is bound to the recoiling SMBH, with 
similar lum inosity as a globular cluster or ultracompa ct dwarf galaxy, and very high velocity 
dispersion (IMerritt et al.ll2009l ; lO'Leary fc Loebll2009l ). The other signa ture is an offset tidal 



flare d ue to the recoiling SMBH tidally disrupting surrounding stars. iKomossa fc Merritt 
( 120081 ) have used both analysis and TV-body simulations to calculate the tidal disruption 
rates of the recoiling SMBH for bound and unbound stars. They found that the tidal disrup- 
tion r ates, in most cases, are smaller than in stationary systems. Recently, lO'Leary fc Loeb 
(120 111 ) have investigated the bound stellar density profile and tidal disruption evolution 
around the recoiling SMBH. They found that the tidal disruption rate will monotonically 
fall as a power-law after the recoil. There is also an anal ytical work focusing on the tidal dis- 
ruption rate during a short period just after coalescence (IStone &: Loebll201ll ). which predicts 
multiple tidal disruption flares in few years or decades after the coalescence. 

Most of the above works about tidal disruption of the recoiling SMBH have focused 
on theoretical analysis and neglected the impact of background stars. Since the real evo- 
lution process in such kind of system is very complex, a detailed study on the recoiling 
SMBH co-evolving with the galactic center is essential. Unlike previous work, we evolve the 
entire system before and after the recoil with unbound stars included. Our work focuses 
on investigating the co-evolution between the recoiling SMBH and surrounding stars with 
tidal disruption and accretion processes. To it come true, we use a special high-accuracy, 
parallel direct iV-body code accelerated by many-core hardware (lp -Grape and <p -GPU) 



( Berczik et al.ll2005l ; lHarfst et al.l 12007; ISpurzem et a 



simplified tidal disruption scheme (IFiestas et al.ll201ll ). Most of our simulations are calcu- 



20091 ; I Just et al.l l201lh . including a 



lated on the laohu GPU cluster in National Astronomical Observatories of China (NAOC). 
With up to million of particles in simulation, we can find out the dynamical evolution of 
stars near the recoiling SMBHs and estimate their tidal disruption rates, which may give us 
some hints for observation. 

Our results show that the accretion of disrupted stars does not significantly impact the 
dynamical evolution of recoiling SMBH. However, it changes the tidal disruption rate com- 
pared to a stationary SMBH system in a galactic center. Most of tidal disruption events are 
contributed by unbound stars when the SMBH passing through the galactic center, and the 
disruption rate in average is ~ 10" 6 M Q yr _1 , which is roughly an order of magnitude lower 
than the stationary case, ~ 10" 5 M Q yr" 1 . Besides, the tidal disruption rate for oscillating 
SMBH far from the galactic center is roughly one or two orders of magnitude lower than 
stationary SMBH. We also find a bound stellar system aro und the recoil i ng: SM BH with 
a mass of ~ 0.7% black hole mass, which is consistent with IMerritt et al.l (120091 ). Except 
for bound stellar systems, we find that there is a "polarization cloud" composed of stars 
dynamically perturbed by the oscillating SMBH. This "polarization cloud" is axisymmetri- 
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cally diffuse in a large spatial area which is beyond the gravitational influence region of the 
oscillating SMBH. For this reason, most of the stars in the cloud are not permanently bound 
to the SMBH. This kind of clouds are the echo of dynamical friction. The calculations for 
parameter dependence show that the variation of initial mass of the recoiling SMBH can 
speed up or slow down the dynamical evolution, and impact the tidal disruption rate. In 
addition, the variation of density slope can obviously change the dynamical evolution and 
tidal disruption rate of the recoiling SMBH. 

We give our galactic model and simulation method in Section [2J Some results about 
bound/unbound stellar systems and tidal disruption rates for the stationary/recoiling SMBH 
are showed in Section [3j The dependence of our results on other parameters are also investi- 
gated. Section H] gives discussion about the observational implications and a short summary. 



GALACTIC MODEL AND iV-BODY SIMULATIONS 



2.1. Galactic Model 



For simplicity, we adopt a spherical Dehnen model to describe the ellipticals or bulges of 



galaxi es ( lDehnenlll993l ). which is different from a Sersic model used in lGualandris fc Merritt 
( 120081 ). In a Dehnen model, the space density profile follows 



p(r) 



7 



Ma 



1) 



4vr rT(r + a) 4 "T 

where a is the scaling radius, M is the total mass of galaxies, and the slope index 7 is a 
constant between interval [0,3). The mass distribution is proportional to r~ 7 for r <C a and 
r -4 for r 3> a. In this model, it is easy for us to deri ve the galactic potential, cumulative 



mass, and half- mass radius respectively (lDehnenlll993l ) 



$(r) 



I ~2^7 

X 



r+a 
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M(r) = M 
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(3) 



ri/2 = a [2 
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where G is the gravitational constant. From Equation fl2]) the escape velocity at radius r for 
7 < 2 is 



I 2$ ( r ) I 

2 GM\ 1/2 
2-7 a ) 



r 



r + a 



2-7 



1/2 



for 7 < 2. 



(5) 



Thus the escape velocity from the galactic center for 7 < 2 is 



Ksc(o) 



2 GM 



1/2 



(6) 



.2-7 a 

If there is a SMBH with mass Mbh at the galactic center, we can estimat e its influence r adius 
r in f by using the definition with stellar mass M*(r ^ r- m {) = 2Mbh as in Merritt ( 120061 )1*1 



a(2M BH /M) — 
1 - (2M BU /M)^ 



(7) 



For simplicity, we adopt model units G = M = a = 1 thereafter. In this new unit, the 
quantities relate to physical quantities with the scale factor of time, velocity, and length, 
respectively 



[T] 



[V] 



[R] 



GM Y 



-1/2 



A/ r \ -1/2 / \ 3/2 

M \ 1 ( 77/2 N 



L « xl ^- 1 » I,1 (i4)" ( 



1 kpc 



yr, 



GM 



1/2 



a J 

655.8 x (2^7 - 1)- 1/2 



t\j \ 1/2 / \ -1/2 



10 11 Mrr 







1 kpc 



kms , 



- = (2^-l)[^) kpc- 



(8) 



(9) 
(10) 



1 This kind of definition is validated only for a singular isothermal sphere nucleus. Here we adopt it also 
for 7 ^ 2 to estimate r- ln { at an order of magnitude. 
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With the model units, Equations fll]), (jEJ), and © can be rewritten 



r 1/2 = [2^) - 



Ksc(o) 



2-7 



1/2 



(12) 



(2M] 



BH 



l 

3-7 



i - (2M BH ; 



(2M] 



BH 



l 

3-7 



3-7 



(13) 



For a fiducial model adopted in the following calculations with 7 = 0.5, M = 4 x 10 M , 
^bh = 4 x 10 7 M Q and ri/ 2 = 1 kpc, we can derive V esc (0) ~ 847 km s^ 1 and r in f ~ 29 pc. 



2.2. Tidal Disruption Scheme in A^-Body Simulation 

Compared to a scattering experiment, direct iV-body simulation can well represent the 
dynamical co-evolution of the recoiling SMBH and surrounding stars, especially for the 
stellar interactions. However, it can not deal with other process well, for example, the tidal 
disruption by SMBH. To solve this problem, we implement a special tidal disruption scheme 
in our simulation. 

A star with mass and radius will be tidally disrupted if it approaches a BH within 
the tidal radius Jmilslll975l : lReeslll988h 



(14) 



Here /1 is a dimensionless parame ter of order u nity. It is roughly 10~ 6 r in f — 10 -5 r in f for the 
real conditions in galactic nuclei ( jMerrittll2006l ). Based on Equation f[T4"j) . to make sure that 
the tidal radius is larger than the event horizon of the BH, there is a limitation for Mbh- 
For a disrupted star with solar radius R and mass M , a Schwarzschild BH should have 
^bh ^ 10 8 M Q . This limitation mass can be heavier when the BH has rapid spin. 

In order to investigate the tidal disruption of the recoiling SMBH with A^-body simu- 
lation, we have to find a proper way to represent the tidal radius in our simulations. As 
shown in Equation (tHj), the tidal radius is proportional to the radius of the disrupted star. 
However, in our A^-body simulations without stellar evolution scheme, stars are point like 
particles, whose radii are not defined. We can use the approximate scaling relation between 
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r m { and r t discussed above to set a proper value for r t , but this value is so small that there 
will be only few tidal disruption events recorded in the simulation. That is because of the 
limitation on particle resolution. Our simulation particle number is not enough even with 
N = 10 6 . 

To solve this problem, a simplified strategy is adopted. We assume a larger r t in our N- 
body simulation to represent the tidal radius. The particle which get close to SMBH within 
r t will be tidally disrupted and removed from the system with its mass added onto the 
SMBH. With this method, we can directly follow the growth of Mbh and calculate the tidal 
disruption rate. However, the large r t we adopted here will overestimate the tidal disruption 
rate. The solution is to carry out a series of simulations with decreasing r t and extrapolate the 
results to the regime corresponding to real galaxies. Our extrapolation scheme is explained 
in Section 13.3.11 in detail. We do not consider the growth of r t for accreting SMBH or the 
momentum transfer from tidal disrupted stars to the recoiling SMBH, because the total 
mass of disrupted stars is negligible compared to the BH. A test simulation taking into 
account momentum transfer shows no significant difference in either the dynamical evolution 
of recoiling SMBH or the stellar disruption rate. 

2.3. Numerical Method 

To investigate the evolution of the recoiling SMBH with tidal disruption carefully, we 
make a series of integrations with different parameters which are listed in Table [TJ The 
first column is the sequence number for different models. Columns (2) - (4) give particle 
number N, tidal radius r t , and initial stellar density slope 7 respectively. Column (5) is the 
initial mass ratio between Mbh and the total mass of system. Column (6) gives initial recoil 
velocity in the unit of escape velocity at the galactic center. 

All of our integrations adopt a parallel direct iV-body <£>-Grape/(/?-GPU code with 
fourth-order Hermite integrator and simplified tidal disruption scheme. They all have time- 
step accuracy parameter 7] = 0.01 and a softening length e = 10 -5 . Most of our models 
adopt iV = 10 6 equal mass particles, with initial mass of SMBH Mbh = 0.001 to represent 
the ellipticals/bulges and central SMBH respectively. To find out the dependency of the 
results on particle number, several sets of integrations with 50K and 250K particles have 
been included. There is also a calculation with Mbh = 0.002 to check the dependency on 
Mbh- As mentioned before, for the purposes of extrapolating to the smaller tidal radius in 
the real galactic center conditions, we have varied the r t from 5 x 10 -5 to 5 x 10~ 2 , which 
is roughly corresponding to ~ 10" 3 r; n f to ~ r- m f for 7 = 0.5. Since the density slope prior 
to the recoil is under debate, we choose 7 = 0.5 as our fiducial value, and also investigate 
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other cases with 7 = 1.0 and 1.5. The recoil velocity is set to Vk = 0.7V esc for most of the 
integrations, and there are calculations without recoiling velocity for com parison. Besides, an 
i ntegr ation with \% = 1.1V^ SC and 7 = 1.0 is included for comparing with lKomossa fc Merritt 
(|2008h . 

Since Dehnen model provides the stellar distribution in a galaxy without central SMBH, 
we have to first to set up consistent initial conditions for our system. At first, we put a SMBH 
into the center of a stellar system with Dehnen profile, and evolve the entire system for 50 
simulation time units to relax the core region. After that the system is dynamically and 
thermally relaxed in relation to the SMBH gravitational potential. At this stage, we do not 
consider stellar tidal disruption process. However, we want to make sure that stars inside the 
tidal disruption loss cone of the central SMBH have already been removed by tidal accretion 
before the recoil of SMBH. Therefore, as a second step, we switch on the tidal disruption 
scheme in the program and run it again with the SMBH still in equilibrium in the center - 
this costs a few dynamical timescales at the influence radius. Thus all particles within the 
tidal disruption loss cone can be removed before we apply a kick to the SMBH. After setting 
up the initial configuration of the stellar distribution of the systems, we artificially give a 
recoil velocity to the SMBH and follow its dynamical evolution. 



3. RESULTS 

3.1. Dynamical Evolution of the recoiling SMBH 

To carefully investigate the dynamical evolution of recoiling SMBH, we have a long time 
integration with model 06 to t = 200, as shown in Figured] It is clear that the evolution of 
the recoiling SMBH can be easily divided into two different phases. The phase I has obviously 
damping trajectory from t = to t ~ 50. After that, the SMBH evolves into phase II, where 
the trajectory damped very slow for a long time. In phase I, the traject ory of recoiling 



SMB H can be well predicted by Chandrasekhar's dynamical friction theory (I Chandr asekhar 



19431 ). However, this stage will just keep till the SMBH's oscillation decayed to the size 
of core, where a slowly decayed phase II begins. That is because the stellar mass interior 
to SMBH's orbit is roughly equal to M B h at the end of phase I, where the assumptions of 
Chandrasekhar's dynamical friction theory are invalid. We have tried to fit the trajectory 



analytically in phase I with the same method as iGualandris fc Merrittl (120081 ) used, and got 
a similar result even though a Dehnen model is adopted here instead of their core-Sersic 
model. 



If we continue the integration, the orbital oscillation of the SMBH will slowly decay and 
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finally achieve to a thermal equilib rium with surrounding stars, which is the so-called phase 
III in iGualandris fc Merrittl (120081 ). Here we have not continue our integrations into that 
stage because phase III is actually similar to the case of a stationary central SMBH, which 
is out of the scope of this paper. We will just focus on phase I and phase II in the following 
integrations. Thus all of our models integrate to t = 100 in order to make sure that the 
evolution can achieve to phase II. 

With reg ard to the orbital mot i on of the SMBH, our results agree qualitatively well 
with those of IGualandris fc Merrittl (120081 ) . That is beca use the growth of the b lack h ole 
due to tidal accretion, which did not take into account in IGualandris &: Merrittl (120081 ) . is 
relatively small. 



3.2. Compact Star Clusters around the Recoiling SMBHs 



A recoiling SMBH can carry a retinue of bound stars. iMerritt et al.l (120091 ) predict that 
there will be a "hypercompact stellar system" (HCSS) bound to the recoiling SMBH. We 
ca n analytically estima te the bound population of the HCSS with Equations (la) and (6a) 



in 



Merritt et al.l (120091 ) . For example, the roughly boundary radius r ej of HCSS should be 



r ej = 



BH 



V 2 



(15) 



0.7Ksc ~ 0.8, we have r e j ps 1.5625 x 10 3 in our simulation 
unit. Thus we can define a factor /b to describe the fraction of bound stars, 



with our Model 06, where V e j 



M 



BH 



*i(7) ( ^ 

\ r mf 



3-7 



(16) 



where Mb is total mass of bound stars. Based on Equations ffl~5]) and (|13|) . we have fb 
3.8845 x 10" 5 Fi(7). For the case of 7 = 0.5, with Fi(0.5) in Figure la of IMerritt et al. 
(120091 ). we have /b ~ 0.001. That means the ratio for HCSS mass to Mbh should be ~ 0.1% 



corresponding to just few particles for 1M simulation. It should be aware that this fraction 
is just a lower limit because our "two times Mbh" assumption actually overestimate r; n f . To 
confirm the conclusion above, we investigate this problem with iV-body simulations. 

In the simulation, it is very difficult to distinguish whether a star is bound to SMBH or 
not. That is because most of star particles with negative total energies relative to the SMBH 
are loosely bound. Only few of them could orbit the SMBH for more than one orbital period. 
For this reason, we identify a star to be bound to the recoiling SMBH if following two criteria 
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are both satisfied: the particle (1) has negative total energy relative to the SMBH and (2) 
remains bound for at least one orbital oscillation period of the recoiling SMBH. We find that 
most of the "bound" particles around SMBH particle will be evaporated before the end of 
simulation. This effect is more significant for our escaping SMBH calculation with model 10, 
where the number of bound particles decreases monotonically. For other oscillating SMBH 
calculations, there are many unbound particles can be captured by the SMBH. That means 
the interaction between the HCSS and stellar background is very important. 

Figure [2] shows the evolution of bound stars around a recoiling SMBH in model 06. 
At t — 0, when V e j = 0.7V esc , there are five bound particles, roughly consistent with the 
estimation from Equation f fl6|) . The number of bound particles, iVbd, slightly increases 
during the evolution from phase I to phase II. Because in phase II (1) the low relative velocity 
between stars and the recoiling SMBH increases the probability of dynamical capture and 
(2) the stellar density around the SMBH in phase II is higher. Our results suggest that a 
HCSS may grow via capturing unbound single stars. However, the low particle resolution 
and large statistical error prevent us from further investigation of the growth of HCSS. 

The HCSS may be detected as an off-nuclear compact system of si milar size and stellar 



mass of a globular clu sters but of having very high internal velocities (IMerritt et al.l 12009 



O'Leary fc Loebl l201ll ). However, that is not the whole story. As shown in the Figures [3] 
and HI we find that there is a stellar cloud composed of stars which are strongly impacted 
by the recoiling SMBH at the first apocenter. The total mass of this cloud is comparable 
to Mbh- They form a quasi- axisymmetric stellar cloud and always have a dense region 
around SMBH. This kind of clouds gives an evidence for the echo of dynamical friction. For 
this reason, we name that stellar cloud as "polarization cloud" . In our simulation of model 
06, we choose a time snapshot when the SMBH arrives at its first apocenter, and pick up 
those polarization cloud members from the E tot — r phase space which includes all of star 
particles. All of polarization cloud members in the E tot — r phase space are outliers around 
SMBH particle and can be easily distinguished. Thus we can collect those polarization 
cloud particles and trace their evolution. However, it is very difficult to find an automatical 
solution to efficiently distinguish those polarization cloud members for every time snapshot. 
To deeply investigate this kind of polarization clouds, an efficient method to extract those 
outlier particles is needed. This problem should be solved in our following work. 

In Figure [31 we plot that polarization cloud at time t = 1.863 Myr, when the recoiling 
SMBH arrives at its first apocenter. All the parameters here are obtained from model 06. 
In order to investigate the origin and evolution of those cloud members, we look back the 
data and follow the evolution of the stars from time t = 0. Figure H] gives the evolution of 
those stars and the recoiling SMBH from t = until t = 42.6 Myr, when our computation 
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stopped. The snapshots with t = Myr, 0.852 Myr, 1.863 Myr in FigureH]show the evolutions 
of the cloud members from beginning toward the first apocenter of the recoiling SMBH. The 
following three snapshots correspond to the time when the recoiling SMBH comes back from 
the first apocenter and passes through the density center for the second time. The snapshot 
at t = 7.984 Myr shows the SMBH passing through galactic center for the third time, and 
the last snapshot is for the end of computation. Our calculation results and Figure [4] imply 
that the shape of polarization cloud varies with the oscillations of the recoiling SMBH. The 
comparison between the first and the last two panels vividly shows the changes of the cloud 
members during the oscillation. 

As shown in Figure [21 there are only small fraction of stars really bound to the SMBH, 
whereas the mass of the polarization cloud is roughly equal to Mbh- Therefore most of 
these cloud members are unbound to the recoiling SMBH. After the recoil of SMBH, those 
polarization cloud members are accelerated and fall behind the oscillating SMBH. During the 
oscillation, some of the stars gradually expand to larger area while others continuously follow 
the SMBH to oscillate around the galactic center, and form an axisymmetrical distribution 
along the velocity direction of the recoiling SMBH. After multiple interactions with oscillating 
SMBH, these stars obtain energy from oscillating SMBH and diffuse to a large area. This 
is a vivid example that an oscillating SMBH transfer its orbital energy to surround stars 
through dynamical friction. 

In order to investigate the observational properties of this polarization cloud, we focus 
on its most compact stage, when the recoiling SMBH arrive at the apocenter for the first 
time. The distribution of those stars are shown in Figure [31 Here we scale the iV-body 
system to our physical fiducial galaxy model with stellar mass M = 4 x 10 10 M and half 
mass radius ri/ 2 = 1 kpc. For such system, the polarization cloud has mass of M csc ~ 10 7 M Q , 
and a size of diameter D ~ 260 pc. Further calculation shows that, its half mass radius is 
~ 83 pc. This size and mass are much larger than a globular cluster and comparable to 
an ultracompact dwarf galaxy. We also calculate the velocity dispersion of this polarization 
cloud. The line-of-sight velocity dispersion along X, Y and Z axes are 185 km s^ 1 , 190kms _1 
and 189kms _1 individually. That is significantly higher than a typical ultracompact dwarf 
galaxy. Comparing with bound HCSS, this polarization cloud has larger size, heavier mass 
and similar velocity dispersion. Since this cloud is mixed with HCSS, it may bring troubles 
to the detection of HCSS. We also estimate the average surface density inside the half mass 
radius of the polarization cloud. The overdensity of the polarization cloud relative to the 
stellar background is only ~ 4%. This may be because of that the apocenter of the recoiling 
SMBH here is not very far from the galactic center, where the surface density of the stellar 
background is still relatively high in model 06 with 7 = 0.5. For a SMBH with a higher 
recoil velocity and a larger apocenter, the situation may be improved. To investigate this 
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problem, more simulations are needed. We will investigate this problem in our future work. 

It should be noted that our example above is just for one time snapshot of the evolving 
SMBH. Every snapshot could have a similar stellar cloud. The composition of these clouds 
may form an interesting structure which can trace the trajectory of the SMBH. The evo- 
lution and observation characters of this kind of structures will be discussed in our future 
work. For a recoiling SMBH with highly enough velocity to escape from the galaxy, that 
kind of polarization cloud may not be distinguishable because many stars will be dropped 
behind the escaping SMBH. Without the multiple interactions between oscillating SMBH 
and surrounding stars, the polarization cloud can not survive for a long time. 



3.3. Tidal Disruption 

3.3.1. Calculate Tidal Disruption Rate Numerically 

For a solar type star disrupted by a SMBH with mass M BH = 4 x 1O 7 M , based on 
Equation (j!4p and fllOp . the tidal radius is ~ 8 x 10~ 6 pc, corresponding to ~ 10~ 8 in simula- 
tion units with our fiducial model. As mentioned in Section I2TTI because of the limitation on 
particle resolution, it is practically impossible with present computing capabilities to calcu- 
late the tidal disruption rate directly (with a realistically small tidal radius) through direct 
iV-body simulations. 

To solve this problem, we perform a parameterized study, which treats both the particle 
number N and the size of the tidal disruption radius r t as free parameters. All of the 
simulations have been done with these varying paramete rs, and physical conclusion s are 



drawn from extrapolating to the real galactic conditions. iFiestas fc Spurzeml (120101 ) have 
shown in Fokker-Planck models that such scaling is reliable and provides physically correct 
results. Besides, tidal disruption events in our simulations, as well as the case in real galactic 
center, do not occur frequently, which means that the increase of Mbh is discrete. Thus it 
is difficult to compute the tidal disruption rate through calculating accreted masses in time 
bins. Instead, we measure the tidal disruption rate in a cumulative, averaged way - fitting 
the mass growth of the black hole and derive its time derivative. 

In order to find out the extrapolation relation, the dependence of tidal disruption rate 
on particle numbers and tidal radius should be investigated. Figure [5] shows the disruption 
rates changes with different N and r t , with parameters are N = 250K, 1M, r t = 5 x 10~ 3 , 1 x 
lO- 3 , 5 x 10~ 4 , 5 x 10~ 5 , 7 = 0.5, M BH = 0.001 and V oi = 0.7Ksc- The left panel shows 
that the dependence of averaged tidal disruption rates on N and r t . Here we calculate the 
averaged disruption rates for phase I, phase II and entire process separately with different 
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N and r t . These integrations show that the dependence of tidal disruption rate on particle 
number is very weak. Further discussions are showed in Section 13.41 



Based on the loss cone feeding theory (IFrank fc Reeslll976l ; iLightman fc Shapiro! 119771 ) 



contrary to empty loss cone case, the tidal disruption rate for full loss cone does not depend 
on N. It should be notice that the tidal disruption rate here is for mass disruption rate, and 
the disrupted particles indeed depend on N. The slight different average rates between 1M 
and 250K cases means that the loss cones in our integrations are nearly full and the changes 
of particle numbers do not bring significant impact. 

Opposed to particle numbers, the changes on r t strongly influence the disruption rates. 
A detailed study is in the right panel of Figure |5l which shows the evolution of disruption 
rate - r t dependence for several special points with integrations for iV = 1M. Five special 
snapshot have been picked out. Here "1st. APO.", "1st. BACK", "1st. DOWN", "2nd 
APO." and "Phase Tran." refer respectively to the snapshot that the oscillating SMBH the 
first time arrives at the apocenter, the first time returns to the density center, the first time 
reaches the opposite apocenter, the second time return to the apocenter and the time of 
transition from phase I to phase II. Our results show that the disruption rates increase for 
every r t value during the evolution of the recoiling SMBH. That is because a decayed orbit of 
SMBH leads to a growing stellar background around it. Besides, the similar linear disruption 
rate - r t relation as left panel has been found. However, the relation here for the first four 
snapshots are not as good as the average results in the left panel. That may because our 
particle resolution for small r t is not good enough in phase I. 

Based on our result data, it seems that there is an approximate linear relation between 
r t and disruption rates. As a result of full loss cone condition, the disruption rates will be 
proportional to cross section S and density p(r) near the recoiling SMBH, 



M ~ 5p(r)U rcla , (17) 

where U re i a is average relative velocity between the recoiling SMBH and the surrounding 
stars. 

It is difficult for us to analytically derive stellar distribution around a recoiling SMBH. 
However, the cross section can be wrote as 

s-«U + ( 18 ) 

which can write the tidal disruption rate as 
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M * nr?p(r)v rela ( 1 + ™%™) . (19) 

V r t v ie\a J 

For r 3> r t that is 2GM B u/r t ^> U^ cla (even if one takes into account that v^ elSL contains a 
contribution from the SMBH kick velocity), so we are in the gravitational focusing regime, 
where 

S K 2 * G J'™ n . (20) 

t; rela 

Thus we have 

27rGM BH 

M ~ — — —p(r)r t . (21) 

f^rela 

If p and w re i a are constant with time, the tidal accretion rate will scale linearly with 
r t . This relation is clearly only correct in a time averaged sense, because when the SMBH 
oscillating in the galaxy, all quantities will be strongly changed within a single orbit (density, 
stellar velocity dispersion and SMBH velocity, the two latter quantities determining u re i a )- 
However, in a running time averaged over several orbits, our results show that there is not a 
large variation in these quantities. So we can define a dimensionless tidal accretion parameter 
a by the following relation: 

tdynj^- = a _ I . (22) 

M BR V re \ a t dyD 

Here tdyn is the dynamic time scale around the position of the recoiling SMBH. Our simu- 
lations imply that across all of our model families a does not depend strongly on Af or r t , 
but only on the evolutionary phase. Thus a can be used in extrapolations to the real system 
with very small r t and large N. Based on the analysis above, though it is complicated to 
get an accurate scaling relation, we still can estimate the tidal disruption rate through linear 
extrapolation of tidal radius. 



3.3.2. Tidal Disruption Rate for the Stationary SMBH 

For comparison purposes, we calculate the tidal disruption rate of a stationary SMBH 
without recoil in the galactic center. A series of simulations with different r t values are carried 
out. Figure [6] shows the static disruption rate for r t = 5 x 10 -3 , 1 x 10 -3 , 5 x 10~ 4 , 5 x 10 -5 , 
corresponding to model 15, 17, 05 and 19 respectively. As mentioned in Section [3.3.11 the 
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tidal radius for our fiducial model is ~ 8 x 10 -6 pc, corresponding to ~ 10 -8 in our simulation 
units. Our simulation results with different particle numbers for stationary case show that 
the tidal disruption rates weakly depend on the particle numbers. Thus we can loosely accept 
a full loss cone condition for the stationary case. As discussed in Section 13.3. 1\ there is a 
roughly linear correlation between tidal disruption rate and rt for full loss cone case. Based 
on Equation (j2T|) . we can extrapolate our simulation results to the real galaxy conditions 
with smaller tidal radii. Figure [6] shows that the tidal disruption rate for r t = 10~ 3 is 
~ 8 x 10~ 6 in our simulation units. That can be extrapolated to r t ~ 10~ 8 , and gives us the 
tidal disruption rate ~ 10~ 10 , which corresponds to the order of magnitude ~ 10~ 5 M Q yr" 1 
for our fiducial model. 



3. 3. 3. Tidal Disruption Rate for the recoiling SMBH 

In a recoiling SMBH system, there are two populations of stars contributed to the tidal 
disruption rate. One is the bound stars around the SMBH, and the other is those unbound 
stars encountered w i th the recoiling SMBH. If the bound stars dominate, as demonstrated in 



Komossa fc Merrittl (120081 ) for an escaping SMBH, the tidal disruption rate should roughly 
keep a constant level at the beginning, and finally drop down because of the consumption 
of bound stars. On the contrary, if the tidal disruption rate is dominated by unbound stars, 
the results will depend on stellar density, relative velocity of SMBH, and the size of tidal 
radius, which described by Equation (TT9|) . 

If we just consider about the contribution of the unbound stars, for an oscillating SMBH, 
its tidal disruption rate near the galactic center should be higher than the apocenter region 
because of the different stellar density. However, w rola can also change the result. For instance, 
we have an estimation in our Model 16 with r t = 5 x 10~ 3 , 7 = 0.5, Mbh = 0.001 and V e j = 
0-7V csc ~ 0.8. Our simulation result shows that the first apocenter of the recoiling SMBH is 
around tbh ~ 1-4. With the assumptions that w rc ia(^ = 1-4) ~ 0.1 and tJ re ia(^ = 0) ~ 0.8, 
we can estimate the disruption rate ratio from the galactic center to the first apocenter 
region. By substituted into Equations (CQ) and (fl9|) . that ratio is ~ 8. When scaling to the 
real galactic center conditions with r t ~ 10 -8 , it should be ~ 3. Moreover, this ratio also 
depends on U re i a . If v re \ a near apocenters can be larger than 0.1, that ratio may be higher. 
Unfortunately, we can not obtain an accurate value for U re i a - For estimation, it can be seen 
that the tidal disruption rates near the galactic center are several times or even an order of 
magnitude higher than the apocenter region. 

That conclusion is confirmed by our simulation results. Figure [7] gives the distribution 
of the tidal disruption events count N t d relative to distance r in model 16. It shows that most 
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of the tidal disruption events appear around the density center, which indicates a boosted 
tidal disruption rate when the oscillating SMBH pass through the galactic center. A SMBH 
around the apocenter with a nearly zero velocity always corresponds to a low u re i a , which 
should give us a high disruption rate. However, as Equation ( TT9]) shown, the low stellar 
density environment makes the tidal disruption events less than at galactic center. The 
higher disruption rate near the galactic center means that, for an oscillating SMBH, most 
of tidal disruption events are contributed by unbound stars. For this reason, unbound stars 
are very important for calculating tidal disruption rate. 

Figure [8] shows the tidal disruption rates for a recoiling SMBH with different r t = 
5 x 10~ 3 , 1 x 10~ 3 , 5 x 10~ 4 , 5 x 10~ 5 , which belong to model 16, 18, 06 and 20 respectively. 
To increase the accuracy of our fitting results for disruption rate in phase I and phase II, 
we neglect the transition boundary between the two phases in the fitting. Therefore, we fit 
disruption rate for phase I and phase II separately. As a result, the fitted results are not 
continuous between two phases. 

The left part of Figure M shows that the tidal disruption rates in phase I are linearly 
increasing. For different r t values, there is also a similar linear relation to stationary SMBH 
case, even though there is an increasing tidal disruption rate instead of a constant value. Here 
the calculation with model 20 (r t = 5 x 10 -5 ) is not included because the particle resolution 
is too low to estimate disruption rate. With the same extrapolating relation described in 
Section 13.3. 1\ we estimate that the averaged tidal disruption rate for phase I should around 
the order of magnitude ~ 1O _6 M yr" 1 , which is about an order of magnitude lower than 
stationary SMBH in the galactic center. 

For the phase II, when the orbit of the recoiling SMBH damped to relative small region, 
the variations of stellar density and velocity of SMBH are smaller than phase I. Thus the 
difference of disruption rate between apocenters and galactic center is not so intensive. Since 
the tidal disruption events during entire phase I are few in our simulation, the accuracy for 
calculating tidal disruption rate is not very good. However, there are plenty of tidal disrup- 
tion events in phase II, which provides a better accuracy to calculate the tidal disruption 
rate in phase II. Actually, as shown in the right part of Figure El there are constant disrup- 
tion rates as similar as the stationary case. That is because (1) the stellar densities around 
SMBHs in phase II and in the stationary SMBH case are similar and (2) their loss cones are 
both full. 
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3.4. Dependence on Particle Number 

Since 1M particles do not have enough high resolution to represent a real galaxy, it is 
important for us to find out the dependence of our results on the particle number. Only with 
this considered, the direct iV-body simulation results for dynamical properties of the recoiling 
SMBHs and surrounding stars can be extrapolated to real conditions in galactic nuclei. To 
bring this forth, a series of integrations with various N and same other parameters have been 
done, which relative to models 01, 02 and 06 with N = 50K, 250K, 1M respectively. Figure [9] 
shows the particle number dependence of the recoiling SMBH orbital oscillation and mass 
increasing for model 01, 02 and 06 respectively. It can be seen that the 250K and 1M cases 
have similar oscillate amplitudes in phase II, whereas the 50K run can not well represent 
the damped evolution of the recoiling SMBH. Besides, only the calculations with 250K and 
1M particles can give relatively smooth mass increasing curve, and they all achieve to the 
similar final mass. 

Limited by the particle resolution, we can not obtain good enough mass evolution data 
to fit the disruption rates in phase I both for 50K and 250K integrations with model 01 and 
02. The only thing can be done is to calculate average tidal disruption rates for phase I and 
II, as discussed in Section 13.31 with Figure [5j However, those problems do not exist for all 
of 1M integrations and other 250K integrations. For this reason, most of our calculations 
adopt 1M particles. As shown in Figure the difference of average disruption rates between 
250K and 1M is very small, which also indicates that the loss cones are nearly full. 

With 1M particles simulation, the mass ratio of SMBH to stars should be ~ 10 3 . This 
ratio is smaller than the real galaxy condition which should be ~ 10 6 — 10 9 . The small 
mass ratio between S MBH and stars may lead to a higher Brownian velocity comparing with 



real galaxy condition (IMerritt et al.ll2007l ). However, most of our simulations are terminated 



before phase III, which makes the amplitude of oscillating SMB H is greater than Brownian 



amplitude. Detailed discussion in iGualandris fc Merrittl (120081 ) with series of calculations 



for different particle numbers shows that the particle mass can not significantly impact the 
dynamical evolution for both phase I and II. 

The limitation of particle number can also influnece the relaxation time T r . For a 
homogenous isotropic d istribution w ith equal-mass stars, the two-body relaxation time in 



our simulation units is (jSpitzerl 119871 ) 



0.34a 3 0.34iVcx 3 , , 

Tr(r) « — - = , (23) 

jO(rJm*lnA p(r) In A 
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where In A is Coulomb logarithm and we can estimate it as 

In A ^ In = In . (24) 

When A" is large enough, for instance A" > 10 6 , the relaxation time will be roughly 
proportional to N/lnN. It means that our A^-body simulation with 1M particles always 
give us a shorter relaxation time comparing with real galaxy condition. This fast relaxation 
in A^-body simulation may bring troubles to our results. If the core region in our simulation 
has relaxed before the recoiling SMBH return, our results will departure from the real galaxy 
condition, because the relaxation time in a real galactic core is always very long. 

To solve this problem, we calculate the relaxation time of the central core region for 
both the simulation model and real galaxy condition. We choose r = r m {, N = 10 6 and 
Af = 4 x 10 10 in model 06 to calculate the relaxation time T r (r; n f). Since we know little 
about the exactly value of velocity dispersion a around r = r; n f, we range the a from 10~ 2 
to 1. For the real galaxy condition, our estimations show that, the T r (r in f) is always longer 
than the half period of the oscillating SMBH. That means the core region with a real galaxy 
condition can not relax before the SMBH return. The same result applied to our simulations 
with a > 0.06 for A" = 10 6 . Based on Equation (Q, 0.06 corresponds to ~ 45kms~ 1 , which is 
smaller than most of the galactic center cases. Thus we can conclude that the core region can 
not relax before the recoiling SMBH return both for simulation and real galaxy condition. 

The limitation of particle number in A^-body simulation may also impacts the results for 
the stationary SMBH. In recoil case, the oscillating SMBH always corresponds to the full loss 
cone status both for the simulation and the real galaxy condition. For a stationary SMBH, 
things will be different. Both the analytical calculations and the numerical simulations for 
long term evolution of tidal disruption fro m stationary SMBH show that the loss cone in 



real galaxy is empty for most of the cases (jFrank fc Reesl 1 19761 ; iLightman fc Shapiro! 11977 



Cohn fc Kulsrudlll978l ; iBaumgardt et all 120041 ) . However, our calculations are just focus on 
a relative short period, which correspond to ~ 10 7 yr for our fiducial model. During this 
short period, the loss cone could be roughly full. This prediction has been confirmed by our 
simulation results. 

Based on discussions above, there are some conclusions: (1) both 250K and 1M integra- 
tions are good enough to represent dynamical evolution of the recoiling SMBH, (2) in order 
to investigate the tidal disruption and co-evolution of the recoiling SMBH with surrounding 
stars, it is better to use integrations with A" = 1M, (3) the small difference on average dis- 
ruption rates between 250K and 1M integrations indicates that the loss cones are nearly full, 
otherwise the results should change obviously with various particle numbers, (4) the particle 
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number we adopt in simulation is less than the real galaxy condition, which can lead to a 
shorter relaxation time scale. However, it can not significantly impact our results. 



3.5. Dependence on other Parameters 



Based on the Chandrasekhar's dynamical friction theory (jChandrasekharl Il943l ). the 
dynamical friction time scale is inversely proportional to the mass of the recoiling SMBH. 
For this reason, a heavier SMBH should has shorter evolution time scale in phase I. In order 
to investigate the dependence of evolution time scale on black hole mass Mbh, we run a 
test model 07 with Mbh = 0.002. Comparing with model 06, we find that model 07 gives a 
shorter oscillation period in phase I, which is roughly two times smaller than that in model 
06. That is consistent with the prediction of dynamical friction theory. Besides, the tidal 
disruption rate in model 07 is roughly two times higher than that in model 06 for both of the 
two phases, which is roughly consistent with Equation (ED). For a galaxy with lighter SMBH, 
which has Mbh < 10 _3 M, there will be a longer dynamical evolution time in phase I and 
smaller tidal disruption rate comparing with our results above. Limited by the calculation 
capability, simulations with the lower mass SMBH will cost so much time that we can not 
investigate them carefully. 

We note that in the simulations with very large tidal radius, for example, r t = 0.05 in 
model 13 and 14, rapid increase on Mbh lead to the declined disruption rates. This significant 
growth of Mbh is artificial, which makes the declination of disruption rate is unphysical. 

When considering different density slopes, the results are quite different. The higher 
of 7 value is, the faster of the system evolved. Our integration with model 12 shows that 
the recoiling SMBH will enter phase II around t — 5, which is more faster than evolution 
of SMBH in model 06 with t ~ 50. Thus our models for 7 = 1.0 and 1.5 actually reach 
phase III at the end of integration. That can be easily understood because a larger 7 means 
a denser cusp in the galactic center. For the same reason, a larger 7 also leads to a higher 
disruption rate. 

Figure [10] shows the disruption rates with different 7 values for both the stationary and 
the recoiling cases. The left panel is for stationary integrations and the right panel is for the 
recoiling SMBH with time duration t = 60 — 100. Here we can not compare the same phase 
for different 7 models because of their different evolving times. The results in the right panel 
are phase II for model 06 and phase III for model 09 and 12. It can be seen that, for both 
panels, there are several times differences between every two adjacent 7 values. 

In order to compare with the work in the literature with higher V e j, a simulation with 
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similar initial parameters comparing with model 06 has been done in model 10, where 7 = 1 
and V^j = 1.1V^ SC ~ 1.56. Our result shows that there are just five particles tidally disrupted 
at the beginning of recoil, and following with zero event recorded for the rest of time. Because 
of the low particle resolution, we can only give a lower limit to the tidal disruption rate. 
For a galaxy with mass M = 4 x 1O 1O M , M BH = 4 x 10 7 M Q and r 1/2 = 1 kpc, the 
kick velocity should be V e j ~ 1000 km s~ 



6 x 10 8 M yr x , which is lower than 



Thus our lower limit for tid al disruption rate is 
1O~ 6 M yr^ 1 in Figure 2 of iKomossa fc Merritt 



(120081 ) with similar conditions. For the bound stars, we find that there are only ~ 10~ 3 Mbh 
stars still boun d to SMBH after the recoil, which is roughly as the same order of magnitude 
as the result in iMerritt et al.l (120091 ) . The small amount of bound stars, low stellar density 
outside and poor particle resolution result in the rare disruption events in computation. To 
solve this problem, a series of simulations with larger particle numbers are needed. 



DISCUSSIONS AND CONCLUSIONS 



Supermassive binary black holes in galactic nuclei may under certain conditions get very 
close and coalesce due to strong emission of gravitational waves. Because of the anisotropy of 
GW emission, the coalescence remnant SMBH experiences a strong recoil velocity depending 
on spins of the original pair of black holes. In this paper we investigate the interaction and 
co-evolution of the recoiling SMBHs with their galactic stellar environments, using very large 
direct iV-body simulations with a simplified tidal disruption scheme. It is the first numerical 
investigation for the evolution of the tidal disruption by an oscillating SMBHs with unbound 
stars included. Moreover, we discover a polarization cloud of stars surrounding the oscillating 
SMBHs, which is not mentioned in the literature. 



As argured by iGualandris fc Merrittl (120081 ). the dynamical evolution of the recoiling 
SMBH can be divided into three phases. During the phase I, the trajectory of the recoiling 
SMB H can be well predicted by Chandrasekhar's dynamical friction theory (1 Chandr asekhar 
19431 ). The phase II begins when the SMBH's oscillation decayed to the size of galactic core, 
and the motion of the SMBH damps slowly. At phase III, the orbital oscillation of SMBH 
slowly decay to achieve a thermal e quilibrium with surrou n ding; stars. In our simulation, 
we obtain similar phase I and II as IGualandris fc Merrittl (120081 ) have got. Besides, our 
results show that the mass fraction of bou nd stars relative to SMBH is ro ughly 0.7% at 



the be ginning. This value is consistent with IGualandris fc Merrittl ( 120081 ) and IMerritt et al. 

hood ). 



In addition to following the dynamical evolution of the recoiling SMBH, we estimate 
the stellar tidal disruption rates by the recoiling SMBHs using direct A^-body simulation 



-22 - 



together with our simple tidal disruption scheme. Our results show that the stellar tidal 
disruptions do not significantly change the dynamical evolution of the recoiling SMBH. The 
stellar disruption rate by the recoiling SMBH in phase II is ~ 10~ 5 M yr^ 1 , which is about 
the same order of magnitude of the stellar disruption rate by the stationary SMBH at a 
spherical galactic center. However, in phase I, the stellar tidal disruption rate is about an 
order of magnitude lower. 

In our simulation, a spherical stellar distribution is assumed. This simplified model 
is not accurate to represent a real galaxy. For most of rea l bulge s/ ellipticals, their stellar 



distributions are non- spherical. As argued by IVicari et al.l ( 120071 ). the orbital decay time 
scales of the recoiling SMBHs are prolonged in triaxial galaxies comparing with spherical 
galaxies. That is because most of the recoiling SMBHs in triaxial galaxies can not return 
directly to their galactic centers, where the influences of dynamical friction are the strongest. 
For the same reason, a non-spher ical dark matter ha lo can also significantly increase the 



decay time of a recoiling SMBH (jGuedes et al.ll2009l ). Since most of the tidal disruption 



events happened during the SMBHs passing through the galactic centers, the triaxiality may 
significantly reduce the tidal disruption rate from unbound stars. To investigate the impacts 
of triaxiality, a good solution is to have series of simulations with triaxial stellar distributions. 
However, it is very complicated to generate a stationary triaxial distribution for iV-body 
simulation, and the prolonged orbital decay time needs more calculation resources. Here we 
just follow previous works to investigate the tidal disruptions from the recoiling SMBHs in 
spherical stellar distribution. Further works with triaxial stellar distribution will be included 
in our next paper. 

Actually, the triaxiality can not change our results significantly because th e influence of 



triax iality is moderate in galaxy models with shallow central density profiles (IVicari et al. 



20071 ). Thus our calculations with 7 = 0.5 will not significantly departure from the triaxial 
model. Besides, in most of our simulations, the maximal displacement of the recoiling 
SMBH relative to the density center is no more than 1 kpc. That keeps the trajectory of 
the SMBH inside the bulge, where the impact of dark matter halo are very weak comparing 
with baryonic matter. For this reason, w e negle ct the influence of the triaxial dark matter 



halo which is mentioned by iGuedes et al.l ( 120091 ) . 



Our simulation for escaping SMBH does not have enough particle resolution to estimate 
disruption rates. That needs much larger iV-body particle number in the future simulations. 
Here our integrations are focus on relatively low recoiling velocity. For the assumption that 
M = 4 x 10 10 M and ri/ 2 = lkpc, recoil velocity in our fiducial model is ~ 600 kms" 1 . 
This velocity produces the oscillating SMBH an amplitude less than the effective radius 
of elliptical/bulge. Our results show that most of tidal disruption events occur when the 
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recoiling SMBH passes through density center, but not exact the original center because of 
the orbital precession. The tidal disruption rates around each apocenter are much lower. 
An estimate for average tidal disruption rate at about apocenters in phase I over several 
periods for model 06 has been done. We calculate the average disruption rate when SMBH 
have large distant relative to density center. The disruption rate at large distance is several 
times or even an order of magnitude lower than the average value over whole phase I, which 
implies that the rate for apocenters is one or two order of magnitude lower than that by 
stationary SMBHs. 



Komossa fc Merrittl (120081 ) argued that off-nuclear X-ray tidal flares could be one of the 
observational signatures for an escaping SMBH. Unfortunately, our simulation of model 10 
with 1M particles and the escapable recoil velocity does not have sufficient resolution to 
estimate the disruption rate. We can just give a lower limit ~ 6 x 10~ 8 M Q yr _1 , which does 
not conflict to their work with similar configurations. Besides, our model 09 with 7 = 1.0 
and recoil velocity ~ 600 km s" 1 corresponds to a tidal disruption ra te ~ 10~ 5 M^ yr~ 



This value is similar to the result of Model 1 with V e j = 500 km s 1 in iKomossa fc Merritt 



( 120081 ). Based on our results, an oscillating SMBH with enough large recoil velocity may also 



contribute to off-nuclear X-ray flares around its apocenters with the low disruption rates. 
An HCSS around the recoiling SMBH in our integrations has similar mass as predicted 



by lMerritt et al.l (120091 ) . However, we can not give more information about the star clusters 
limited by particle resolution. In addition, we find a larger cloud of stars around the recoiling 
SMBHs, which are not permanently bound to the black hole but keep a kinematic correlation 
to it for more than an orbital time. It can be described as polarization cloud reacting to 
the gravitational disturbance by the recoiling SMBH, in other words, the echo of dynamical 
friction. The size and mass of the polarization clouds are comparable to ultracompact 
dwarf galaxies, except significantly higher velocity dispersion. The existence of such kind of 
system may bring difficulties for distinguishing the HCSS. More detailed properties of the 
polarization cloud will be studied in our following work. 

Our simulations for different particle numbers show that both 250K and 1M integrations 
are sufficient in following the dynamical evolution of the recoiling SMBH, but the former do 
not have enough resolution for investigating the tidal disruption and co-evolution of SMBH 
with surrounding stars. Different initial masses of the recoiling SMBHs can change the 
results moderately. We find that a higher galactic density slope can significantly speed up 
the dynamical evolution of the system, and the tidal disruption rate is also boosted. 
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Fig. 1. — Trajectory of the recoiling SMBH with model 06. Here the bottom horizontal axis 
and the left vertical axis give the evolved time and the displacement of BH rt,h respectively, 
with the units G = M = a = 1. The right and the top axis are transformed to real 
physical units with the assumption that M = 4 x 10 10 M Q and r\i<i = 1 kpc. It can be easily 



disting uished for the phase I and II which are similar to the result in iGualandris &: Merritt 

(boosh . 
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20 40 60 80 100 

t [Dehnen] 

Fig. 2. — The evolution of bound star numbers around a recoiling SMBH based on the model 
06. The horizontal and vertical axis are the same as Figure [31 the right horizontal axis is 
the number count of the bound particles. The grey dotted line gives the trajectory of the 
recoiling SMBH, and the black solid line gives the variation of the bound particles. 
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Fig. 3. — Distribution of stars which are just strongly impacted by the recoiling SMBH 
at its first apocenter with model 06 for time snapshot t = 1.863 Myr. We project the 
stellar distribution onto X — Z plane, with the physical units under the assumption that 
M = 4 x 10 10 M Q and ri/ 2 = 1 kpc. The large black spot marks the position of the recoiling 
SMBH. It can be seen that all of the impacted star particles form a cloud structure around 
the SMBH at that moment. 
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Fig. 4. — Evolution of stars which are just strongly impacted by the recoiling SMBH at its 
first apocenter with model 06. The units and definition of axis are the same as Figure [3J 
The black circle and the grey dots mark the recoiling SMBH and the stars respectively. The 
black arrow represents the size and orientation of BH's velocity. It can be seen that many 
unbound stars are influenced by the SMBH. 
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Fig. 5. — Dependence of the tidal disruption rate on N and r t . The left panel shows the 
dependence of the average tidal disruption rates on N and r t . The solid line and dashed 
line are the average disruption rates for 1M and 250K integrations respectively. The square, 
triangle, and circles respect to the average disruption rates for phase I, phase II and the 
entire process respectively. The integrations are based on model 06 with different particle 
numbers and tidal radius. The integration for r t = 5 x 10~ 5 with iV = 250K is not included 
because the amount of resulted tidal disruption events is too small to calculate an average 
rate. The right panel shows the disruption rate - r t dependence for several key points with 
integrations for N = 1M. Here "1st. APO.", "1st. BACK", "1st. DOWN", "2nd APO." 
and "Phase Tran." are corresponding to the moment which SMBH the first time arrive at its 
apocenter, the first time come back to the density center, the first time arrive at apocenter 
on the other side, the second time arrive at its apocenter and the transition from phase I 
to phase II respectively. For both of two panels, horizontal and vertical axis give the tidal 
radius and tidal disruption rate in the units ofG = M = a = l respectively. 
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Fig. 6. — Tidal disruption rates for the stationary SMBH with different r t . Here the hor- 
izontal and vertical axis give the evolved time and tidal disruption rate in the units of 
G — M — a — 1 respectively. The dotted, dashed, solid and dash dotted lines are calcula- 
tions with different tidal radii r t = 5 x 10~ 3 , 1 x 10 -3 , 5 x f 0~ 4 , 5 x f0~ 5 , which correspond 
to model 15, 17, 05 and 19 respectively. 
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Fig. 7. — The distribution of the cumulative tidal disruption events count N td relative to r 
with model 16. Here the horizontal and vertical axis give the cumulative tidal disruption 
events count and r respectively. The shaded histogram marks the cumulative tidal disruption 
events count for the different distance. 
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Fig. 8. — The time evolution of the tidal disruption rates with different r t . The left part 
and right part are the fitted results for phase I and phase II respectively. The shaded 
region corresponds to a transition area from phase I to phase II. Axis are labeled as same as 
Figure [6j The dotted, dashed, solid and dash dotted lines are the calculations with different 
tidal radii r t = 5 x 10~ 3 , 1 x 10~ 3 , 5 x 1CT 4 , 5 x 10~ 5 , corresponding to model 16, 18, 06 and 20 
respectively. Since the model 20 with r t = 5 x 10~ 5 does not have enough particle resolution 
to calculate the disruption rate during phase I, it is absent in the left part. 
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Fig. 9. — The particle number dependence for the recoiling SMBH orbital oscillation and 
the mass increasing with model 01, 02 and 06. Axis are labeled as same as in Figure [TJ The 
dotted, dashed and solid lines are the calculations with particle numbers N = 50K, 250K, 1M 
respectively. Left panel is the damped orbital evolutions of recoiling SMBHs with various 
N. It can be seen that the 250K and 1M cases have a similar oscillate amplitude in phase 
II, while the 50K is different from others because of the poor particle resolution. The right 
panel shows the mass increasing for the recoiling SMBH. Here the vertical axis represents to 
the mass of the recoiling SMBH. The different initial mass of the recoiling SMBH is due to 
the special relaxation scheme described in Section 12.31 Only the calculations with 250K and 
1M particles can give us relatively smooth mass increasing curve, and all the calculations 
achieve to similar order of magnitude for final mass. 
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Fig. 10. — The tidal disruption rates for different 7 values and different models. The left 
panel is the disruption rates of stationary SMBHs with model 05, 08 and 11. The right 
panel shows disruption rates of the recoiling SMBHs for the time interval t = 60 — 100, in 
model 06, 09 and 12. Here the horizontal and vertical axis stand for the evolving time and 
the disruption rates respectively, with the units of G = M = a = 1. The solid, dashed and 
dotted lines stand for the integrations with 7 = 0.5, 1.0, 1.5 respectively. 
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Table 1. Parameters of simulation models 



Model No. 


A ! 


rt 


7 


Mbh/M 


Vej / Vcsc 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


01 


50K 


5 x 10~ 4 


0.5 


0.001 


0.7 


02 


250K 


5 x 10~ 4 


0.5 


0.001 


0.7 


03 


250K 


1 x 10~ 3 


0.5 


0.001 


0.7 


04 


250K 


5 x 10~ 3 


0.5 


0.001 


0.7 


05 


1M 


5 x 10~ 4 


0.5 


0.001 


0.0 


06 


1M 


5 x 10~ 4 


0.5 


0.001 


0.7 


07 


1M 


5 x lO" 4 


0.5 


0.002 


0.7 


08 


1M 


5 x 10~ 4 


1.0 


0.001 


0.0 


09 


1M 


5 x 10~ 4 


1.0 


0.001 


0.7 


10 


1M 


5 x 10~ 4 


1.0 


0.001 


1.1 


11 


1M 


5 x 10~ 4 


1.5 


0.001 


0.0 


12 


1M 


5 x 10~ 4 


1.5 


0.001 


0.7 


13 


1M 


5 x 10~ 2 


0.5 


0.001 


0.0 


14 


1M 


5 x 10" 2 


0.5 


0.001 


0.7 


15 


1M 


5 x 10- 3 


0.5 


0.001 


0.0 


16 


1M 


5 x 10" 3 


0.5 


0.001 


0.7 


17 


1M 


1 x 10- 3 


0.5 


0.001 


0.0 


18 


1M 


1 x 10-3 


0.5 


0.001 


0.7 


19 


1M 


5 x 10~ s 


0.5 


0.001 


0.0 


20 


1M 


5 x 10~ s 


0.5 


0.001 


0.7 



Note. — Col.(l): Model sequence number. Col. (2): Particle 
numbers adopted in calculations. Col. (3): Tidal disruption radius 
rt in the units G = M = a = 1. Col. (4): Density slope 7. Col. (5): 
Initial BH mass in the unit of total mass. Col. (6): Initial recoil 
velocity in the unit of escape velocity. All the simulations set e = 
10~ 5 and n = 0.01. 



